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Abstract. We model the dynamics of magnetization in an artificial analog of 
spin ice specializing to the case of a honeycomb network of connected magnetic 
nanowires. The inherently dissipative dynamics is mediated by the emission, 
propagation and absorption of domain walls in the links of the lattice. These 
domain walls carry two natural units of magnetic charge, whereas sites of the 
lattice contain a unit magnetic charge. Magnetostatic Coulomb forces between 
these charges play a major role in the physics of the system, as does quenched 
disorder caused by imperfections of the lattice. We identify and describe different 
regimes of magnetization reversal in an applied magnetic field determined by the 
orientation of the applied field with respect to the initial magnetization. One of 
the regimes is characterized by magnetic avalanches with a 1/n distribution of 
lengths. 



1. Introduction 

Spin ice [Tl[2] is a frustrated ferromagnet -with Ising spins tliat possesses rather peculiar 
properties. First, as a consequence of strong frustration, it has a massively degenerate 
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ground state and retains a finite entropy density even at very low temperatures [3]. 
Second, its low-energy excitations are neither individual flipped spins, nor domain 
walls, but are point defects acting as sources and sinks of magnetic field H [HIS]. The 
concept of magnetic charges, while not exactly new [HIITIIH], has proven very useful in 
elucidating the static and dynamic properties of spin ice [9l [TOl [iTJ [121 US] • It is worth 
noting that these objects are magnetic analogs of excitations with fractional electric 
charge found in the familiar water ice |14j . 

Artificial spin ice is an array of nanomagnets with similarly frustrated interactions. 
The original system made by Schiffer's group had disconnected elongated islands (80 
nm by 220 nm laterally and 25-nm thick) made of permalloy and arranged as links 
of a square lattice [15]. Later versions included a connected honeycomb network of 
flat magnetic wires (THl HZl HSl US] , in which the centers of the wires form a kagome 
lattice, hence the sometimes used name "kagome spin ice" [T^. Whereas it had been 
originally intended as a large-scale replica of natural spin ice, it became clear very 
soon that artificial spin ice has a number of its own peculiar features. For example, 
because the magnetic moments in artificial spin ice are extremely large, on the order of 
10* Bohr magnetons, the energy scale of shape anisotropy due to dipolar interactions, 
10^ K in temperature units [20], effectively freezes out thermal fiuctuations of the 
macrospins meaning that the system is not in thermal equilibrium. Dynamics of 
magnetization has to be induced by the application of an external magnetic field |15) . 
Elaborate experimental protocols involving a magnetic field of varying magnitude and 
direction [21 have been proposed to simulate thermal agitation invoking parallels with 
fluidized granular matter. It remains to be seen whether the induced dynamics yields 
a thermal ensemble with an effective temperature. The analogy with granular matter 
is further reinforced by recent observations of magnetic avalanches in the process of 
magnetization reversal [TH [22] . 

In this paper we present a model of magnetization dynamics in artificial spin ice 
subject to an external magnetic field. Two sets of physical variables are used: an Ising 
variable cr = ±1 encodes the magnetic state of a spin, whereas an integer q quantifies 
the magnetic charge of a node at the junction of several spins. Magnetization dynamics 
are mediated by the emission of domain walls carrying two units of magnetic charge 
from a lattice node, their subsequent propagation through a magnetic element, and 
absorption at the next node. We specialize to the case of kagome spin ice, in which 
magnetic elements form a connected honeycomb lattice [T8j [El [ITl [19] . The model 
can be readily extended to other geometries and lattices with disconnected magnetic 
elements [T31 [531 [21] ■ Some of the results presented here have been outlined 
previously [25] . 

2. Basic features of the model 

Our model is specialized toward an experimental realization described previously |17j . 
That artificial spin ice is a connected honeycomb network of permalloy nanowires with 
saturation magnetization M = 8.6 x 10^ A/m and the following typical dimensions: 
length I = 500 nm, width w — 110 nm, and thickness t = 23 nm. Three nanowires 
come together at a vertex in the bulk. At the edge of the lattice, a vertex may have 
one or two links coming in. 
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Figure 1. Magnetization variables aij = ±1 (arrows) live on links ij of the 
honeycomb lattice. Charges qi = ±1, ±3 live on nodes i. 



2.1. Basic variables: magnetization and magnetic charge 

We label nodes of the lattice by a single index i and nanowires connecting adjacent 
nodes by the indices of its two nodes, ij. In equilibrium, the vector of magnetization 
M points parallel to the long axis of the wire, so we can encode the two states of 
a nanowire by using an Ising variable aij = ±1. In our convention, aij = +1 when 
the vector of magnetization points from node i into node j. This definition implies 
antisymmetry under index exchange, aij ~ ^'^ji- 

We define the dimensionless magnetic charge at node i as 

qi^^<^3i, (1) 
i 

where the sum is taken over the three neighboring sites j. This definition is quite 
natural: since magnetic induction B = /xo(H + M) is divergence-free, the magnetic 
charge Qi of node i equals the flux of magnetic field H out of the node, which in turn 
equals the flux of magnetization M into it: 

Q, = ^ H • dA = - ^ M • dA = -Mtw ^ a,j = Mtwq,. (2) 

j 

Thus Qi is indeed magnetic charge measured in units of Mtw. 

The Bcrnal-Fowler ice rule [2 enforcing minimization of the absolute value of 
charge \Qi\ is usually justified from the energy perspective: the magnetostatic energy 
of spin ice can be written as the energy of Coulomb interaction of magnetic charges. 

The dominant second term — the charging energy of a node — forces minimization of 
magnetic charges in natural spin ice. The "capacitance" C is determined by the 
dipolar and exchange couplings energies of adjacent spins [5]. 
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Although we will see below that these energy considerations are not relevant to 
artificial spin ice within our model, for the moment we will simply adopt the result 
to it. In honeycomb ice, where the coordination number is 3, dimensionless charge qi 
can take on values ±1 and ±3. Minimization of node self-energy would select states 
with 

q^ = ±1. (4) 

Indeed, triple magnetic charges have never been observed in our samples of honeycomb 
ice. Ladak et al. [18l[T9] have found nodes with triple charges. The difference is likely 
due to a higher amount of quenched disorder arising from random imperfections of 
the lattice in the samples of Ladak et al. 

We will find it convenient to use the following notation. A site with a unit charge 
Qi = ±1 has two majority links with aji — qi and one minority link with aji — —qi- 
For site i in Figure [1] the minority link is ij. 



2.2. Basic dynamics: emission of a domain wall 

To reverse the magnetization in a nanowire, one must apply a sufficiently strong 
external magnetic field. The reversal begins when one of the nodes, say i, emits 
a domain wall (w) into link ij, Figure [^a). If the link initially has magnetization 
fTij = ±1, a domain wall can traverse it from i to j only if it has charge of the right 
sign, i.e., — '2.(Jij — ±2. Once the domain wall passes through the link, Uij changes 
its sign. Now a domain wall with the same charge can only traverse the link in the 
opposite direction. 

The critical field He, at which a domain wall is emitted from a node, can be 
estimated as follows [27]. Suppose a node with magnetic charge = ±1 emits a 
domain wall with magnetic charge g^ = ±2 [SI Conservation of magnetic charge 
means that the charge of the site turns to g^ = Tl- The emission process can thus be 
viewed as pulling a charge g^ = ±2 away from a charge of the opposite sign g^ = =pl. 
The maximum force between the two charges is achieved when the separation between 
them is of the order of their sizes a, which is roughly equal to the width of the wire w: 
Pmax = Mo|QiQw|/(47ra^). This force must be overcome by the Zeeman force applied 
to the domain wall by the external magnetic field, i^ext = Mo | Qw | -ffext . Hence the 
estimate of the critical field, 

^ ^ \Q± ^ Mtw ^ Mt_ 
Atto'^ 47ra^ At:w 

For the system parameters used in our previous work Il7, and listed above, this 
estimate yields /lo-ffc = 18 mT. The critical value observed experimentally [28] is 
35 mT. 

One can envision another possible process, wherein the reversal is triggered when 
a site with charge g^ = ±1 emits a domain wall of charge gw — t2 and change its 
charge to qi — ±3. Considerations along the same lines as above show that the critical 
field required to pull apart charges qi = ±3 and gw — t2 is 3i?c- As we will see below, 
magnetization reversal in samples with low quenched disorder occurs well before the 
external field has a chance to reach this value. This explains why triple charges are 
never generated as a result of the emission of a domain wall. 

The estimate for the critical field was obtained under the assumption that the 
external magnetic field Hcxt is applied along the link into which the domain wall is 
emitted. When the field makes angle 6 with the link, it is reasonable to suppose that 
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Figure 2. Magnetization reversal in a single link. At the end of the reversal, the 
domain wall encounters a node with magnetic charge of the opposite sign (a) or 
of the same sign (b). In panel (b), the emission of the domain wall from the left 
node and its propagation along the horizontal link are omitted for brevity. 

only the longitudinal component of the field i?ext cos 9 pulls the domain wall away 
from the node. We thus expect the following angular dependence of the critical field: 

HciO) ^ He/ cose. (6) 

As we will see later in Sec. |31 our educated guess is almost right and that Eq. © 
requires only a minor correction: the angle 9 should be measured not from the axis of 
the link but from a slightly offset direction. This effect is caused by an asymmetric 
distribution of magnetization around a node, which was missed by the simplified, 
mesoscopic model of this section. 

2.3. Basic physics: absorption of a domain wall 

Once a domain wall is emitted into link ij, it quickly propagates to the other end of 
the link, toward node j. Theoretical and experimental studies of domain wall motion 
in permalloy nanowires [IHl [SD] show that walls move at speeds of the order of 100 
m/s in an applied field of just 1 mT. This corresponds to a propagation time of the 
order of 10 ns, which is too short to be observed in most current experimental setups. 
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When the domain wall reaches the opposite end of the link ij, its further fate 
depends on whether the magnetic charge at node j has the same or opposite sign of 
magnetic charge. We consider the two cases in turn. 

If the domain wall and node j at which it arrives have opposite charges, 
q„ = ±2 = —2qj, as in Figure [2ja), the domain wall is attracted to the node. It 
is easily absorbed by the node, whose charge changes to qj = ±1. A new domain 
wall with the same charge gw — ±2 may be subsequently emitted into one of the 
adjacent links jk if two conditions are met: (i) the link has the right direction of its 
magnetization, — '^o'jk and (ii) the external field is sufficiently strong to trigger the 
emission. 

Note that condition (ii) is sensitive to the orientation of the field relative to link 
jk. It also rests on an implicit assumption that the critical field for a new domain wall 
is not affected by the just completed absorption of the previous one. This assumption 
is reasonable if the dynamics of domain walls are strongly dissipative and the energy 
generated during the absorption process is quickly dissipated as heat. Experiments 
with domain walls in nanowires indicate that they possess non- negligible inertia [8], 
and therefore our assumption of strongly overdamped dynamics may not be fully 
justified. Nonetheless, for the sake of simplicity, we shall assume that the dynamics 
are strongly dissipative and that the extra energy brought by the arrival of a domain 
wall does not by itself cause the emission of another domain wall from the same node. 

Consider now the other case, where the domain wall and the arrival node have 
charges of the same sign, q^ = ±2 — 2qj, as in Figure [2ljb). The two charges now 
repel and the repulsion grows stronger as the domain wall approaches the node. 
Under the assumption of overdamped dynamics, the wall stops when the Coulomb 
repulsion between the charges reaches the level of the Zeeman force from the external 
field. One might think that this may be an equilibrium situation, but we show 
as follows that this is not the case. The arriving domain wall generates a strong 
field at the node, whose magnitude is easy to estimate. Since the domain wall is in 
equilibrium, the force applied to it by the external field, F = /zq | Qw I -f^c , is balanced 
by the Coulomb repulsion of the node. By Newton's third law, the domain wall 
applies an equal force to the node. The field created by the wall at the node is 
H = F/\fioQj\ = \Qw/Qj\Hc = 2Hc. This field is added to the externally applied 
field He- The resulting field is sufficiently strong to trigger the emission of another 
domain wall from the node. (This works for any relevant direction of the applied field.) 
The charge of node j changes sign, qj = =f1 = —q„/2, and subsequently absorbs the 
stopped domain wall. 

2.4- Basic physics: quenched disorder 

Imperfections of magnetic links and junctions create local variations of the critical field 
He- If the variations of He result from a large number of small errors, one expects a 
Gaussian distribution of critical fields p{Hc) with a mean He and a width SH^ given 



In the limit of strong disorder, when the distribution width 5 He is comparable to the 
average J?c, nodes with the highest critical fields may fail to follow the scenario shown 
in Figure [2ljb) and remain in a state with a triple charge until the field becomes strong 



by 




(7) 
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enough. Nodes with triple charges have been observed by Ladak et al. [TF1[TH]. In 
contrast, other samples have never shown triply charged defects |17) . indicating that 
these samples are in the low-disorder limit, 5Ht. 0.04//c }28| . 

The distribution width 5Hc can be compared to another characteristic field 
strength, the magnetic field generated by an adjacent node, Hq — Mtw / [AttI'^) . With 
the aid of Equation ([5]), we estimate 

i/o/^c = {a/If « {w/lf. (8) 

If Hq <^ SHc, the Coulomb fields produced by adjacent and more distant nodes can 
be ignored to a first approximation. The Coulomb contribution to the net field on a 
given site is small, but occasionally the redistribution of magnetic charges on nearby 
sites may trigger the emission of a domain wall if the net field is close to the critical 
value. See Sec. 14. II for further details. In the opposite limit, Hq >• SHc, these internal 
fields must be taken into account. The reversal of magnetization on one link alters the 
magnetic charges on its ends. The resulting increments of the total magnetic field at 
nearby nodes, of order Hq, may be sufficient to trigger the emission of domain walls 
from them. Samples we studied previously jl71 128) appear to be in the regime where 
Hq and SHc are comparable. 

3. Microscopic basis for the model 

To test the basic model of magnetization dynamics presented in Sec. [21 we performed 
numerical simulations of magnetization dynamics in a small portion of the honeycomb 
network by using the micromagnetic simulator OOMMF (31j . 

The typical numerical experiment involved a junction of three permalloy magnetic 
wires of length I = 500 nm, width w = 110 nm, and thickness t = 23 nm [17]. We used 
the two-dimensional version of the oommf code with cells 2 nm x 2 nm x 23 nm. (The 
lateral size of the unit cell should not exceed the minimal length in the micromagnetic 
problem, the magnetic exchange length obtained from exchange and dipolar couplings. 
In permalloy, it is about 5 nm [29].) The magnetization field M(r) was allowed to 
relax to an equilibrium state with magnetic charge q = ±1 at the junction, Figure [3] 
An external magnetic field was then applied in a fixed direction and its magnitude was 
slowly increased keeping the system in a state of local equilibrium. Eventually, the 
magnet reached a point of instability when a domain wall was emitted from either the 
central junction or one of the peripheral ends of the wires, depending on the direction 
of the applied field. The wall then propagated to the opposite end of the link reversing 
the link's magnetization. Using those orientations of the field for which a domain wall 
is emitted from the junction, we determined the dependence of the critical field H on 
the angle 9 between the field and the link in which the reversal occurs. Figure lU 

Two features of the angular dependence in Figure S] stand out. First, H{d) is not 
an even function of the angle 9, and contrary to our expectations, the critical field is 
not at its lowest when the field is parallel to the link. Second, the critical fields for 
three different links in the experiment have the same shape but differ in the overall 
scale He- 

We have traced the physical origin for the asymmetric dependence of the critical 
field H{6) to an asymmetric distribution of magnetization at the junction. Figure [3] 
The energetics of the emission process shown in the figure can be described in the 
language of collective coordinates [32] . The soft mode associated with the emission of 
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Figure 3. Reversal of magnetization in a magnet consisting of three joint links in 
an applied magnetic field (vertical arrow). In panels (a) through (c), the strength 
of the field slowly increases from to a critical value as the magnetization adjusts 
adiabatically. In panels (c) through (f), a domain wall detaches from the node and 
quickly propagates through the vertical link; the field value remains essentially 
unchanged. Micromagnetic simulation (oommf). 
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Figure 4. The dependence of the critical field H on the angle 9 between the 
field and the link. The lines are best fits to Equation II12I I. Links 1, 2, and 3 
had He = 53.6 mT, 54.7 mT, and 55.3 mT and a = 19.3°, 19.4°, and 19.4°, 
respectively. The same numerical experiments were repeated three times, with 
the field initially lined up with Link 1, 2, or 3, and then rotated through 180° +6 
from that direction. 
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a domain wall into the vertical link is the domain wall displacement X along the link. 
To the first order in the applied field H and to the second order in X, the energy is 

UiH, X) = C/(0, 0) - fXoXiQ^^H, + Q^yHy) + kX^ /2, (9) 

where Qxxi Qxy a-nd k are phenomenological constants. Generally speaking, the off- 
diagonal component Q^y does not vanish unless the magnetization distribution is 
symmetric under the reflection y i— )• — y. The equilibrium position of the wall depends 
on the direction of the applied field H — {H cosd, H smd,0) as follows: 

^cq = {l^a/k){QxxHx + QxyHy) = {no/k)QHcos{d - a), (10) 

where the offset angle a and effective charge Q are defined through 

Qxx^Qcosa, Qxy^Qsina. (11) 

According to Equation ([TUl) . the relevant component of the magnetic field H is found 
by projecting the field onto the easy axis of a (majority) link, which is rotated 
through angle a toward the minority link. These considerations suggest the following 
modification for the postulated field dependence of the critical field ([5]): 

Hc{e) = HJ cos{e - a). (12) 

As FigureHlshows, this equation provides a good description of the angular dependence 
of the critical field with the offset angle a ~ 19°. The overall scale of the critical 
field He showed variations reffecting small imperfections of links in the simulation. 
For instance, the square lattice of magnetic moments used in oommf simulations is 
incommensurate with links pointing at 60° to a lattice axis and creates edge roughness. 
This observation confirms the proposed model of disorder introduced in Sec. 12.41 

4. Numerical simulations 

The heuristic considerations of Section [2] and the micromagnetic simulations of 
Section [3] suggested a coarse-grained model of magnetization dynamics in which the 
basic degrees of freedom are Ising variables of magnetization (T,j on links and magnetic 
charges qi on sites of the honeycomb lattice. Each link has its own fixed critical field 
He- The critical fields form a Gaussian distribution Q of width 5Hc around the 
mean H^. The average, = 50 mT, was chosen on the basis of our micromagnetic 
simulations, whereas the relative width was set to SHc/Hc = 0.05, a value inspired 
by our experimental observations |28| . Simulations were performed in a rectangular 
sample with 937 links. The edge consisted of "dangling" links with no other links 
attached to their external ends. We choose the initial state with a maximum total 
magnetization that can be obtained by placing the system in a strong magnetic field 
along one set of links. Figure [SKa). Simulation details are described in [Appendix A 

Following initialization, the external field is switched off and reapplied along a 
different direction, at an angle 9 to its initial orientation. Figure [SKb-d). To stimulate 
magnetization dynamics, the rotation angle must be large enough so that H would 
have a negative projection onto at least some of the majority links. When \0\ is 
between 30° and 90°, only one of the three sublattices of links will reverse. Two 
sublattices reverse when 16*1 is between 90° and 150°. The entire lattice undergoes a 
reversal when 16*1 > 150°. 

Aside from the number of active sublattices, there are marked differences in the 
dynamics of the reversal. For small angles of rotation, \6\ < 131°, the reversals 
occur in a gradual and uncorrelated manner, with individual links switching when the 
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Figure 5. Magnetization reversal in an applied magnetic field, (a) The system is 
initially magnetized in a strong horizontal field, (b-d) The field is then switched 
off and applied at 120° to the original direction with a gradually increasing 
magnitude. Open arrows denote links with reversed magnetization. 

applied field reaches the link's critical field. For larger angles, \9\ > 131°, we observed 
avalanches in which long chains of links reverse magnetization simultaneously. 
This kind of switching happens when the sublattice whose magnetization is most 
antiparallel to the applied field cannot switch first because it consists entirely of 
minority links and must wait for one of the other sublattices to begin its reversal. If 
that happens in a higher field, the former sublattice acts like a loaded spring, making 
the reversal nearly instantaneous. A diagram depicting different regimes as a function 
of the field rotation angle 9 is shown in Figure [S] 

In the simplest case, the reversal of magnetization in a link occurs when the 
magnetic field reaches the critical value for that link. The links thus reverse on an 
individual basis, largely independently from the others (but see below). To be more 
precise, the two ends of a link have different critical fields and the reversal begins 
from the end with the lower critical field and stops at the other end. The effective 
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Figure 6. Left panel: Regimes of magnetization reversal. < < 30°: no 
reversal. 30° < 6» < 90°: sublattice B only. 90° < < 131°: B, then A. 131° < 
d < 150°: A and B reverse together. 150° < < 180°: A and B, then C. Similar 
regimes obtain for negative 9, with sublattices B and C exchanged. Right panel: 
Illustration to Equation {TsJ. The Gaussian distribution exp (—x'^/2)/\/2tt (blue 
squares), the modified distribution erfc(a;/\/2) exp (-x^/2)/V2^ (red circles), 
and the best Gaussian approximation, exp {—{x — S)-^ /2iy^) / \/2na with the mean 
S = -0.54 and width a = 0.82 (solid line). 



probability density of the critical fields thus changes from a Gaussian distribution to 

/•oo 

f'{H,) = 2p{H,) j AHp{H) 

' "xpr-<5^%,.fr^?^v ,13) 



It can be seen in Figure [S] (right panel) that the resulting distribution is very close to 
a Gaussian with renormalized mean and width, 

li^^H^- O.MSHc, SH'^ = 0.82(5iJc. (14) 

In our simulation, the renormalized values are H'^ = 48.7 mT and 6H'^/H'^ = 0.042. 



4.1. 30° <9 < 131°.- gradual reversal 

With the field rotated through 6 = 120°, two sets of links have a negative projection 
of magnetization onto the field. In Figure [5{b), they are the horizontal minority 
links and the majority links parallel to the field. Because emission of a domain wall 
into a minority link requires a very high field, it is the majority links that undergo 
magnetization reversal first. The field makes an angle a ~ 19° with their easy axes, 
so the reversal is expected to occur around the field Hi — H'^/ cos (—19°) — 51.5 mT. 

Magnetization reversal in the links parallel to the field alters the magnetic charges 
on all sites. Figure [5Uc). As a result of this change, horizontal links join the majority 
and become capable of reversing their magnetization. The external field makes 
an angle 60° — a « 41° with their easy axes, so their magnetization reversal is 
expected to occur when the field reaches a higher value, H2 = H'^^/ cos 41° = 64.5 
mT. In the presence of disorder, the reversal regions are expected to have finite 
widths, SHi/Hi = 5H2/H2 = 5H[./H[.. For a Gaussian distribution of critical fields, 
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Figure 7. Magnetization reversal curve M{H) in an applied field rotated 
through 120°. Left: Simulated magnetization curve M{H) (red circles) is well 
approximated by the theoretical curve I I15II (solid black line). Inset: semi-log plot 
of the number of avalanches as a function of their length. Right: Experimental 
magnetization curve M{H) (red circles) |28) and the best fit to Eq. Ijl5|l (solid 
black line). 



magnetization measured along the applied field is expected to be a superposition of 
error functions: 

M{H) 1 /H-Hi\ 1 JH-H2\ 1 , , 

— — = -erf ^ + -erf ^ + -. 15 

M^ax 2 \SH1V2J 4 \SH2V2J 4 

The three terms reflect the contributions of the three groups of links with different 
orientations. 

The simulated dependence M{H) is shown in Figure [7] along with the theoretical 
curve (|15|) that takes into account the renormalization of the Gaussian distribution 
parameters p^ . 

A close inspection of the simulated curve M{H) shows that on occasion several 
adjacent links reverse simultaneously due to a positive feedback during the reversal. 
When magnetization of a link is reversed, magnetic charges at its ends are switched. 
The net magnetic field on an adjacent site, projected onto its easy axis, increases by 

Ai/ = 2iJocos41° - (2i7o/3)cosll° = 0.86iJo = 0.74 mT. (16) 

The extra field is not negligible on the scale of the critical-field distribution width 
5H'^ = 2.0 mT. It can help to stimulate the emission of a domain wall at an adjacent 
site if that site's critical field is not too high. This kind of positive feedback causes 
avalanches, in which magnetization reversals occur nearly simultaneously in links 
residing along a one-dimensional path determined by the orientations of easy axes. 
For example, an avalanche occurring in the background of a fully magnetized state of 
Figure [SJb) would travel along the vertical direction. In the limit of small feedback, 
/S.H <C 5H'^, the distribution of avalanche lengths is exponential. Indeed, if the link 
starting an avalanche of length n has a critical field i7, n — 1 of its neighbors must 
have critical fields in the range between H and H + IS.H . The probability to find such 
a collection of links is 

P„ ^ n 1 [p[H)^Hr-^p{H) dH = " ' (17) 
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for a Gaussian distribution of critical fields ([7]). The distribution of avalanches seen in 
the simulation is shown in the inset of Figure [7] along with the theoretical distribution 

m- 

These results can be directly compared to the experimental reversal curve 
measured in the same geometry [28], Figure [7] (right panel). Although the overall 
scale of the magnetic field is substantially lower, the data are well fit by Eq. ([T5l) with 
Hi = 35.9 mT and H2 = 45.9. The ratio of the reversal fields, H2/H1 = 1.28, agrees 
well with the theoretical value H2/H1 = cos (— 19°)/ cos41° — 1.25. The relative 
widths are dHi/Hi = 0.037 and 6H2/H2 = 0.046. 

The magnetization curve M{H) was also measured experimentally [28] and 
simulated for = 100°, with similar results. The experimentally measured reversal 
fields were Hi = 34.7 mT and H2 = 91.5 mT and relative widths SH^/Hi = 0.033 
and SH2/H2 — 0.047. The reversal field ratio was H2/H1 — 2.64 in the experiment, 
somewhat off the theoretical value H2/H1 = cos 1°/ cos 61° = 2.06. 

Overall, it appears that our model provides a reasonably good description of 
magnetization reversal when the field is reapplied a.t 9 — 120° to the direction of initial 
magnetization. In this regime, the reversal proceeds in two well-defined stages, each 
involving one subset of links. During each stage, links reverse largely independently, 
although sometimes the reversal in one link changes the field on a nearby site and 
triggers magnetization reversal there. The reversal fields are given approximately by 
the equations 

Hi^ II'J cos {120° ~9~ a), H2 ^ H'J cos {180° - 6 - a) . (18) 

The reversal follows the two-stage scenario as long as Hi < H2, oi 9 < 150° — a = 131°. 
For larger field rotation angle 9, the reversal proceeds in a very different manner. 

4-2. 131° < 9 < 180°.' reversal with avalanches 

When the field is rotated through 9 — 170° relative to the direction of magnetization, 
the theory described in Section l4Tl no longer applies. Because H2, the reversal field 
of horizontal links, is lower than Hi, these links should reverse first. However, that 
is impossible because in the initial, fully magnetized state. Figure [8](a), these are 
minority links whose critical field is roughly 3Hc (Section 12.21) . i.e., much higher than 
H2 = i?c/cosll° « He- For this reason, a horizontal link does not reverse until 
one of its neighbors, a majority link, reverses and in the process alters the charge at 
one of the horizontal link's ends. This converts the horizontal link into a majority 
link enabling it to reverse magnetization. It turns out that this mode of reversal is 
accompanied by long magnetic avalanches. 

In the simplest scenario, the dynamics begins with the reversal of the weakest 
link with the critical field near Hi, Figure |8Kb). The reversal turns the horizontal link 
next to it into a majority link, which is now ready to reverse since the applied field 
exceeds its critical field: H Hi > H2. A q — —2 domain wall emitted from its left 
end travels to the right end where it encounters a site with charge —1, Figure [8](b). 
As discussed in Section 12. 3[ the arriving domain wall induces the emission of another 
domain wall into an adjacent link. Figure [2{b). The magnetization of that link gets 
reversed, bringing us to the state shown in Figure |8l[d) . The cycle repeats creating an 
avalanche. In effect, we have a q — +2 charge moving along a zigzag path parallel to 
the applied field and reversing magnetization of the links along the way. The process 
continues until the moving charge reaches the edge of the system so that an avalanche 
extends from edge to edge. 




Figure 9. Simulated magnetization reversal curve M{H) in an applied field 
rotated through 170°. Inset: a log-log plot of the number of avalanches versus 
their length (red circles) and a fit to the power-law distribution Ijl9|l (solid line). 
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A different scenario may take place if the system has "weak" hnks that trigger 
the reversal when the applied field is at or below H2- These can be links at the edge 
of the system or some sort of defects. Their reversal converts one of the horizontal 
links (critical field Hd) to the majority status, as shown in Figure E^b). When the 
applied field reaches a value sufficient to induce the reversal of that link, an adjacent 
link is also reversed as described above, Figure [DJc-d) . The next horizontal link down 
the line (critical field Hc2) will switch immediately if Hd > Hc2- The switching will 
continue until the avalanche comes to a stubborn link whose critical field exceeds Hd- 
Its reversal will happen in a higher applied field, possibly triggering another avalanche. 

If the first reversal occurs in a link whose critical field Hd is at the lower end 
of the critical field distribution, the first avalanche will be short because it is unlikely 
that a large number of subsequent links will have even lower critical fields. As further 
avalanches get terminated at links with higher critical fields, their lengths will tend 
to increase. Toward the end of the reversal, avalanches will begin with links whose 
critical fields are near the higher end of the distribution. These avalanches will be 
particularly long. The last avalanche in a given string of links will terminate at the 
edge or will meet an avalanche traveling in the opposite direction. These qualitative 
considerations anticipate a wide distribution of avalanche lengths. Indeed, we show in 
[Appendix B| that the avalanches have a power- law distribution of lengths. 



Remarkably, this result applies to any distribution of critical fields, not just a Gaussian 
one, and numerical simulations confirm this picture. 

As can be seen in Figure El magnetization reversal begins in an applied field 
H H2 - 6H2 ^ 47 niT, where H2 is given by Eq. ([HI). At that point, the 
reversals include single pairs of links from two sublattices. Long avalanches, involving 
as many as n = 10 and more links, are observed by the time the applied field reaches 
H K, H2 + SH2 = 51 mT. The length distribution is well fit by the power law ([T9l) as 
can be seen in the inset of Figure [HI 

The third sublattice reverses in much higher fields, H w Ht, = 77 mT, where 



This stage of the reversal proceeds in the gradual manner described previously. 
5. Discussion 

The dynamics of magnetization in artificial spin ice is a complex problem. In this 
paper, we have presented a simple model for this system in terms of coarse-grained 
physical variables (Figure [T]), Ising spins aij living on the links of the spin-ice lattice 
and magnetic charges qi residing on its sites. Inspired by our earlier studies of magnetic 
nanowires [331 132) , where magnetization reversal is mediated by the propagation of 
domain walls, we have expressed the magnetization dynamics in spin ice in similar 
terms. Magnetization reversal in individual links of the lattice proceeds through 
the emission, propagation, and absorption of domain walls with magnetic charge 
= ±2. Coulomb-like interactions between the magnetic charges of the walls 
and lattice sites play a major role in the dynamics. For example, the magnitude 
of the critical field, required for the emission of a domain wall, is set by the strength 
of magnetostatic attraction between a domain wall and the magnetic charge of the 
lattice site. These heuristic considerations have been confirmed and refined through 



Pn - C/n. 



(19) 



H3 ^ H'J cos {2A0° ~ e - a), 



(20) 
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micromagnetic simulations of a small portion of the spin-ice lattice containing a few 
links. 

Quenched disorder is another major element affecting the magnetization 
dynamics. Small imperfections of the artificial lattice are expected to produce a 
Gaussian distribution of critical fields. The experimentally measured curve [IQ is 
consistent with a Gaussian shape and width SHc/Hc « 0.05. 

The dynamics of magnetization reversal strongly depends on the direction of 
the external magnetic field. If the field is applied at a small angle relative to the 
magnetization of a (fully magnetized) sample, 9 < 131° for the parameters we used, 
the reversal proceeds in a gradual way, with links reversing more or less independently 
of each other, when the strength of the applied field exceeds the threshold of a given 
link. For larger angles of rotation, the reversal proceeds in one-dimensional avalanches 
that can easily span the entire length of the system. The reversal in one link with a 
critical field H triggers the reversal in several others along the chain. The avalanche 
stops when it encounters a link whose critical field exceeds H. In this regime, avalanche 
lengths are distributed as a power law, P„ = C/n. 

It should be pointed out that we model the magnetization dynamics in artificial 
spin ice as a purely dissipative process, in which the system moves strictly downhill 
in the energy landscape. Such a picture is very different from an earlier approach 
extending the notion of an effective temperature to these far-from-equilibrium systems 
[20l |34] . Whereas energy of a microstate plays a major role in the effective thermal 
approach, our method puts the focus on energy gradients, or forces between magnetic 
charges. 

This study has a limited scope. We focus on a continuously-connected honeycomb 
network realized in several experimental studies [HI [HI [H] and cover only the basic 
regimes of its magnetization dynamics. Figure [6l Interesting phenomena arise at 
the boundaries between different regimes, particularly when the field is completely 
reversed, 9 = 180°. In this case, avalanches lose their unidirectional character and 
become random walks. As the magnetization reversal proceeds, avalanches can begin 
to intersect and block one another. 

Our method can be easily extended to connected networks with other geometries 
such as square spin ice [15]. Budrikis, Politi and Stamps used a similar heuristic 
approach to study the dynamics of disconnected magnetic islands |35) . 
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Appendix A. Simulation procedure 

For a given applied external field, the total magnetic field H for each site is computed 
as a sum of the applied field and the Coulomb fields generated by the charges at the 
neighboring sites and domain walls (see Section . For simplicity, we only include 
the fields from first and second-neighbor sites. Fields of further neighbors decrease 
rapidly and tend to oscillate in sign. For each link attached to a given site, the 
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program checks whether the net field has a negative projection — H cos {9 — a) 
onto the Hnk's easy axis, Eq. p^ . If < 0, the program calculates the weakness 
of the site and link, W = \Hc\ — He- The site and link with the largest W in the 
sample are considered to be the weakest. As the applied field increases, the largest 
W becomes positive, triggering the emission of a domain wall from the weakest site 
into the weakest link. The domain wall propagates to the other end of the link where 
it is absorbed, either immediately or after the emission of another domain wall as 
described in Section [2] Once the reversal process that started with the weakest site is 
complete, the program looks for the next weakest site. The process is repeated until 
there are no positive W in the system. Spin ice rules are satisfied at each site at all 
times. No thermal effects are considered. 



Appendix B. Statistics of avalanches in the presence of a weak link 



Here we derive the statistics of avalanches discussed in Section 14.21 In this case, 
the reversal begins in Link 1 (critical field H) and spreads to consecutive Links 
2,3, ... ,n (of the same sublattice) as long as their critical fields are lower than H. 
The avalanche stops when it encounters link n + 1 whose critical field exceeds H. 
The probability density of the critical-field distribution is p{H) and the cumulative 
probability distribution is 

P{H) = p{H')dH'. (B.l) 



Consider an avalanche beginning on link k with a critical field between H and 
H + dH. The k — 1 preceding links must have critical fields less than H. If the 
avalanche has length n then links k + l,k + 2,...,k + n— 1 must have critical fields 
less than H, whereas link k -\- n must have a higher critical field. The probability of 
such a distribution is 

f^{H)dH = [P{H)f-^ p{H)dH[P{H)r-^[l- P{H)]. (B.2) 

However, if the avalanche terminates on link L, the last link of the chain, the factor 
1 — P{H) drops out because there is no link L + 1: 

ff^^^H) dH = [P(i?)]^-" p{H) dH [P{H)r-\ (B.3) 

The probability to find an avalanche of length n is found by summing this distribution 
over the initial position of the avalanche k and integrating over the critical field H . 
Performing the sum first, wc find 

L — n-t-l L — n 

/„ ^ ^ fn= PP^^^ + p{P^+''^^ - = (B.4) 

k=l k=l 

The integration of the resulting expression yields the expected number of avalanches 
of length n, 

/oo 
fn{H) dH 
-OO 

/oo pi 
[P{H)Y'-^p{H)dH ^ / P"-MF = l/n. (B.5) 

Note that Fn is an expectation number of avalanches, not a probability distribution 
normalized to 1. Observing an avalanche of length n does not exclude the possibility 
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of observing an avalanche of a different length n' in the same chain during the same 
reversal process. The probability that an avalanche will have length n is 

L 

Pn=Fn/Y,Fn^l/{n\nL) (B.6) 

n=l 

for large L. 
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